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In evolutionary dynamics, the probability that a mutation spreads through the whole population, 
having arisen in a single individual, is known as the fixation probability. In general, it is not possible 
to find the fixation probability analytically given the mutant's fitness and the topological constraints 
that govern the spread of the mutation, so one resorts to simulations instead. Depending on the 
topology in use, a great number of evolutionary steps may be needed in each of the simulation 
events, particularly in those that end with the population containing mutants only. We introduce 
two techniques to accelerate the determination of the fixation probability. The first one skips all 
evolutionary steps in which the number of mutants does not change and thereby reduces the number 
of steps per simulation event considerably. This technique is computationally advantageous for some 
of the so-called layered networks. The second technique, which is not restricted to layered networks, 
consists of aborting any simulation event in which the number of mutants has grown beyond a 
certain threshold value, and counting that event as having led to a total spread of the mutation. 
For large populations, and regardless of the network's topology, we demonstrate, both analytically 
and by means of simulations, that using a threshold of about 100 mutants leads to an estimate of 
the fixation probability that deviates in no significant way from that obtained from the full-fledged 
simulations. We have observed speedups of two orders of magnitude for layered networks with 10 000 
nodes. 

PACS numbers: 87.23.Kg, 89.75.Fb, 02.10.Ox, 02.50.-r 
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I. INTRODUCTION 

We consider directed networks in which each node is in- 
habited by a single individual of a population and whose 
edges represent the possibilities for an individual's off- 
spring to replace some other individual. Such networks 
provide the substrate on which the evolution of the pop- 
ulation can be studied given the constraints imposed by 
their structures. In this modality of evolutionary dynam- 
ics, known as evolutionary graph dynamics since its in- 
troduction in [1], the population evolves in discrete time 
steps, each of which involves the fitness-based selection of 
an individual for reproduction and the use of its offspring 
to replace one of its out-neighbors in the network. The 
chief quantity one targets in such studies is the proba- 
bility that a mutation arising at a randomly chosen in- 
dividual, henceforth called a mutant, of the otherwise 
homogeneous population eventually spreads through all 
the population. This probability is known as the fixation 
probability. 

In the last decade, the study of several other phenom- 
ena has been approached from a similar perspective of 
interacting agents. Such phenomena have included dif- 
ferently constrained forms of the dynamics of evolution 
[2, 3], the spread of epidemics through populations [4], 
the emergence of cooperation in biological and social sys- 
tems [5|-[8|, and various others |9l-[ll|. In most cases, what 



the interacting agents do, driven by either competition or 
the goal of promoting cooperation, is to spread informa- 
tion through the network in order to attempt to influence 
the states of other agents. In general, network struc- 
ture is a major player in affecting the global outcome 
of such interactions, and this holds to the extent that 
subtle structural changes can have relevant consequences 
[ll, [l2l - [l7| . The importance of network structure, in fact, 
is also central in several other areas, as for example those 
discussed in [18-20]. 

The fixation probability is very heavily influenced by 
the structure of the underlying network as well. In rare 
cases it is possible to calculate it analytically from both 
structure and the relative fitnesses of the individuals 
[l|, l2ll - [2^ , but in general one has to resort to simulations 
of the evolutionary steps. Given the network and the 
mutant's fitness, the simulation is conducted as a num- 
ber of independent events, each of which starts by placing 
the mutant at a randomly chosen node and then carries 
out the evolutionary steps until either fixation (all nodes 
contain mutants) or extinction (no node contains a mu- 
tant) occurs. The fraction of events ending in fixation is 
an estimate of the fixation probability. This simulation- 
based approach to obtaining the fixation probability can 
be very time-consuming, not only because many indepen- 
dent events are needed, but also because each event can 
require a significantly large number of steps to converge. 
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FIG. 1: (Color online) Probability densities associated with 
the number of steps required for extinction and fixation of 
a mutant 10% fitter than the remainder of the population. 
Data are log-binned to the base 1.5 and represent averages 
over 10^ events for the X-funnel, 10^ events for the random 
networks (10^ events for each of 10^ graphs; in the case of the 
random network with a Poisson-distributed number of out- 
neighbors, each of these graphs has at least 95% of the nodes 
in the GSCC). Arrows indicate the means. All the networks 
have 1 555 nodes. 



This is illustrated in Fig. [TJ which suggests two proper- 
ties of the simulation process. The first is that events 
ending in extinction usually require substantially fewer 
steps to conclude than those ending in fixation. The 
second is that the number of steps required for fixation 
varies widely with network topology. 

The various network topologies used in Fig. [1] recur 
throughout the paper, so we pause momentarily to intro- 
duce them. They are the iC- funnel [1], a random gen- 
eralization thereof ^], and the directed variant of the 
Erdos-Renyi random graphs [26| discussed in [23, l28| . 
For b^K > 1 integers, the K-funne\ of base b has node set 
partitioned into K subsets, called layers and numbered 
0,l,...,i^— 1, such that there are b^ nodes in layer k. 
The if-funnel, therefore, has {b^ — l)/{b— 1) nodes. For 
k = 1,2,. ..,i^ — 1, an edge exists directed from each 
node in layer k to each node in layer k — 1. Additionally, 
an edge exists directed from the single node in layer 
to each of the nodes in layer K—1. The i^- funnel is 
generalized by any network in which nodes occupy layers 
and edges exist only from all nodes in a layer to all nodes 
in the next according to some cyclic arrangement of the 
layers. The generalization we use [25] is constructed ran- 
domly by first placing one node at each layer and then 
adding one extra node at a time to a layer that is chosen 
randomly with a bias that favors those layers that are 
nearest to the layer having the most nodes "upstream" 
from them. The other random networks we use [23, [28| 
are such that a randomly chosen node has a number of 
out-neighbors given by the Poisson distribution of mean 
z > 1 [so a giant strongly connected component (GSCC) 
is expected to exist with high probability]. The simula- 



tion of the evolutionary dynamics on such a network is 
confined to its GSCC. 

We also remark, before proceeding, that some of our 
own recent work related to the fixation probability has 
been strongly influenced by the computational difficulties 
associated with estimating it. For example, in [25| we set 
out to create a randomized growing mechanism for lay- 
ered networks that would result in fitness amplifiers in 
the sense of [1]. The mechanism we created gives rise 
to the i^-layer random networks introduced above and, 
indeed, we were able to demonstrate the desired ampli- 
fication effect for 1 000 nodes. However, we fell short 
of demonstrating to full completion that the same holds 
for significantly larger networks, due mainly to the very 
large times required to compute the fixation probabilities 
of the thousands of candidate networks. 

Our central concern in this paper is the devising of 
simulation strategies that can make the calculation of 
the fixation probability substantially faster while main- 
taining accuracy or reducing it only imperceptibly. We 
proceed along two tracks. The first is targeted at elim- 
inating the evolutionary steps at which no change is ef- 
fected to the number of mutants in the network. This 
occurs whenever the simulation prescribes that a mu- 
tant's offspring is to replace another mutant or that a 
non- mutant's offspring is to replace another non- mutant. 
Our results in this track are better suited to the layered 
networks introduced above; we demonstrate them for the 
K-funnel. The second track we pursue builds on the re- 
alization that, if in general it takes a lot more steps for 
fixation to occur than for extinction, then detecting early 
in the course of a simulation event that fixation is highly 
likely to occur can be used as a surrogate to the eventual 
detection of fixation and thereby reduce the number of 
necessary steps. We have found that, nearly regardless 
of the network for a large number of nodes, there exists 
a threshold number of mutants beyond which fixation is 
practically guaranteed to take place. We give results for 
a wide variety of networks. 

We organize the remainder of the paper in the following 
manner. In Section [IT] we briefly review the key notions 
related to the fixation probability. Then in Sections IIIII 
and II VI we pursue the two tracks outlined above for com- 
puting the fixation probability more efficiently, respec- 
tively for layered networks and for unrestricted networks. 
We conclude in Section IVl 



II. FIXATION PROBABILITY 

Let P be a population of n individuals and, for i e P, 
let Pi C P\{i} he the set of individuals that an offspring 
of i can replace during the evolution of P. Let also Qi 
be the set of individuals whose offspring can replace i. 
These give rise to a directed graph, called D, whose set 
of nodes is P and whose set of edges, denoted by E, 
contains the edge (i^j) if and only if j G Pi (equivalently, 
i G Qj). Each individual i has a fitness fi>0 associated 



with it, and similarly to each edge (i, j) there corresponds 
a probability Wij such that ^^^p. Wij = 1 for all i G 
P. The dynamics of evolution that we consider occurs 
in a sequence of steps. At each step, an individual i is 
chosen with probability proportional to /i, then another 
individual j G Pi is chosen with probability Wij^ and 
finally j is replaced by an offspring of i having fitness fi . 

The fixation probability of D^ denoted by p, is the 
probability that a mutation spreads through all of P 
given that it arises at one single individual and that all 
individuals in the remainder of the population have the 
same fitness. The value of p depends on the structure 
of D and on the ratio r of the mutant's fitness to that 
of the other individuals. Moreover, it is the relationship 
between p and r that determines whether evolution is 
driven primarily by natural selection or by random drift: 
essentially, natural selection predominates when p and 
r are highly correlated, random drift otherwise. Note 
in this context that, if D is not strongly connected (i.e., 
there exist nodes i and j such that no directed path leads 
from i to j), then p > if and only if there exists an in- 
dividual from which all others are reachable. This may 
cause random drift to be the main driver of evolution, so 
henceforth we assume that D is strongly connected (thus 
p > necessarily). 

This type of evolutionary dynamics can be described 
by a discrete-time Markov chain of states 0, 1, . . . , n, each 
representing a possible number of mutants in D. In this 
chain, states and n are absorbing and all others are 
transient. If 5 is a transient state, then from s it is possi- 
ble to move to state s -\- 1 (with probability Ps), to state 
s — 1 (with probability Qs), or to remain at state s (with 
probability 1 — Ps — Qs)- If we denote by Pn-\-i{n \ s) the 
probability that, having started at state 5, this (n + 1)- 
state system eventually enters state n, then it is well- 
known that 



Pn+i{n I s) 



1 + En^l nr=l ^v/Pv 



(1) 



(cf. [29| and references therein), whence p, which is given 
by Pn+i(^ I 1), is such that 



1 



^^T.u=lllv=lQv/Pv 



(2) 



When the probabilities Wij are such that ^a^q. 



w 



J* 



1 for alH G P (i.e., not only do the probabilities associ- 
ated with the outgoing edges of i sum up to 1, but also 
those of the incoming edges), and only then, the isother- 
mal theorem of [l], [30| establishes that Pv/Qv = ^ for all 
i; G {1, 2, . . . , n — 1}. In this case, and assuming that 
r 7^ 1 (i.e., the mutation is either advantageous or disad- 
vantageous, but not neutral), it follows from Eq. (J2j) that 
the fixation probability, now denoted by pi, is 



Pi 



1-1/r 
1 - 1/r^' 



(3) 



This includes the Moran process [2], in which D is such 
that Pi = P\ {i} for alii e P and Wij = l/\Pi\ for 
all (i, j) G E (we use \X\ to denote the cardinality of 
set X). It also includes the more general case in which 
the Wij^s are thus constrained and, moreover, all nodes 
have as many incoming edges as outgoing edges, provided 
this number is the same for all nodes. Other noteworthy 
cases that also employ these locally uniform Wij^s are 
the directed graphs that in [1] are shown to be fitness 
amplifiers with respect to the expression for pi. One 
example is the if-funnel, for which substituting r^ for r 
in Eq. (|3]) yields p = Pk as n ^ oo, where 



PK 



1 - 1/r^ 
1 - l/r^n 



(4) 



III. FIXATION IN LAYERED NETWORKS 

We henceforth use Wij = l/|Pi| for all (i, j) G E exclu- 
sively. For t > an integer, let E{t) C E be the set of 
edges (i, j) such that, at the tth evolutionary step, node 
z is a mutant but node j is not. If P^{t) denotes the 
probability that, at step t + 1, the number of mutants in 
D increases (necessarily by 1), then 



p^it)= E 



{iJ)^E{t) 



n + M(t)(r-l)' 



(5) 



where M{t) is the number of mutants at step t [so M(0) = 
1]. Letting mi{t) be an indicator of whether node i is a 
mutant at step t [i.e., mi{t) = 1 in the affirmative case, 
otherwise], this expression can be rewritten as a sum 
over all edges in E: 



p+it) = 



itS 



m,it)[l-mj{t)] 



n + M(t)(r-l) ^^ ^ 
r 



n^M(i){r- 1) 
mi{t) 



E 



E 



mi{t)mj{t) 



■ (6) 



And since D is strongly connected by assumption, we 
obtain 



P+(t) 



n + M(t)(r- 1) 



M{t)- Y. 

iiJ)eE 



mi{t)mj{t) 



(7) 

Similarly, the probability P (t) that the number of mu- 
tants decreases (necessarily by 1) at step t + 1 is 



p-{t) 



1 



n + M(t)(r-l) 
mj{t) 



E 

ihj)eE 



E 

iid)eE 



mi{t)mj{t) 



■ (8) 



The number of mutants remains unchanged at step t-\-l 
with probabihty 1 - P+(t) - P~{t). 

These expressions for P~^{t) and P~{t) ahow the 
evolutionary dynamics to be observed from some in- 
teresting perspectives. For example, any D for 



which 



is such that 



. m 



it)m\ 



p = pi IS sucn mat I](j jjgB . 
E,eP"^iWE,gQ, 1/l^d = Ej^P^jit) = M{t), and 
therefore P~^{t)/P~{t) = r = Ps/Qs for any state s G 
{1, 2, . . . , n — 1} such that M{t) = s. This, essentially, 
is the argument behind the isothermal theorem. More- 
over, for r sufficiently close to 1, what makes P^it) and 
P~{t) differ from each other is the balance between M{t) 



and J2(i j)eE'^ji'^)/\-^'i\^ where the former does not de- 
pend on the topology of D (given t) while the latter does. 
Thus, for example, if we consider the i^- funnel and let 
Mk{t) be the number of mutants in layer k at step t, we 
have 



E 

(i,J)eE 



mj{t) 



2^ lyk 



fe=0 

b[M{t)-MK-i{t)] 



MK-l{t) 



.(9) 



Readily, maintaining a relatively high value for the ratio 
P^ {t) / P~ {t) (a strong forward bias) in the case of the K- 
funnel depends crucially on how close M{t) and Mx-iit) 
are to each other, i.e., on how close the largest layer is 
to containing a significant fraction of the M(t) mutants. 
The fact that P^{t)/P~{t) ^ r^ as n ^ oo indicates 
that the topology of the i^-funnel is in fact successful at 
maintaining the necessary distribution of mutants. 

Equations (O and (J8]) are also useful in that they pro- 
vide an alternative mechanism for simulating the evolu- 
tionary dynamics. Instead of repeatedly choosing z, then 
j G Pi to receive z's offspring, until the mutation either 
spreads through the whole of D or dies out, we use the 
two equations to decide, at each step t, whether the num- 
ber of mutants will increase, decrease, or remain the same 
at step t + 1. In the former two cases we choose the nodes 
to be involved, create or destroy a mutant as the case may 
be, then compute P+(t + 1) and P~(t + 1). By doing so, 
all steps in which no mutant is created or destroyed are 
skipped. Of course, in order for this alternative to be 
computationally attractive the suppression of these steps 
has to compensate for the additional effort to calculate 
the probabilities at every step that is not suppressed. 

While it is unclear that this will be so in the general 
case, for layered networks like the i^-funnel, in which all 
nodes in the same layer are topologically identical to one 
another (they all have the same in- and out-neighbors), 
the simulation can be conducted by keeping track only 
of the number of mutants in each layer and the alterna- 
tive becomes attractive. We then consider the generaliza- 
tion of the iC- funnel obtained by letting layer k have any 
number n/c > of nodes. In this case, we can rewrite 
Eqs. ([7]) and ([8j) by decomposing each of P^it) and 
P~{t) into K summands, each related to a pair of subse- 
quent layers. Thus, we obtain P^{t) = Xl/c^ -^k(^) ^^^ 
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where Mk{t) is as for the K- funnel. In these expressions, 
and henceforth, layer numbers are incremented or decre- 
mented modulo K. 

For layered networks such as these, once it has been 
decided that a mutant is to be created or destroyed at 
step t + 1, the layer at which this is to happen is k with 
probability proportional to P^{t)^ in the case of creation, 
or P^it)^ in the case of destruction [consistently with 
this, notice that P+(t) = if Mk+i(t) = or Mk{t) = 
Uk, and that P~ {t) = if Mk{t) = or Mk+i{t) = 
n/e+i]. And once Mk{t-\-l) has been updated from Mk{t), 
it follows from Eqs. ([lO]) and ([11]) that only P+(t + 1) 
and Pfe+_i(t + 1), or p-{t + 1) and Pj;_^{t + 1), need to 
be calculated. This is so because, although M{t -\- 1) is 
updated from M{t) as well, the only effect this has is 
to alter the factor in the denominator of Eqs. (p!Q|) and 
(pTj) that is common to all layers. As a consequence, the 
probabilities corresponding to any layer other than k or 
k — 1 need not be calculated. 

Computational results on the i^-funnel are given in 
Fig. [2] for a variety of n and K values, where we show 
the speedup afforded by the use of Eqs. (p!Q|) and (pTj) to 
compute the fixation probability. This speedup is defined 
in terms of processor time and indicates, in all cases ex- 
amined, a reduction to less than half the time required 
by the simulation that goes through all the evolutionary 
steps. We note, however, that the approach we discuss 
in Section HVl allows for much more significant speedups. 



IV. FIXATION IN ARBITRARY NETWORKS 

In the absence of the computational facihtation pro- 
vided by layered networks, which as we have seen allows 
the fixation probability to be computed more efficiently 
by skipping all steps of the simulation in which no mu- 
tant is created or destroyed, for an unrestricted topol- 
ogy we turn to the alternative strategy of attempting an 
early stop of each simulation event based on how many 
mutants there are. The central question is whether there 
exists a threshold number of mutants which, once crossed 
from below, ensures that fixation is bound to occur with 
probability as close to 1 as one wishes. We provide an 
affirmative answer in what follows. 

The probability that the mutation eventually dies out, 
given that s mutants are originally present, equals 1 — 
P^+i(n I s). We denote it by Qn+i(0 I -5), and it follows 
from Eq. (jT]) that 

}^u=s llv^l Qv/Pv 



10 



Qn+l(0 I S) 



^^Eu^lUv=lQv/Pv 



(12) 



We are interested in the probability that, conditioned on 
the fact that extinction does actually occur, the number 
of mutants eventually grows from the initial s to some 
fixed value M G {s,s + l,...,n — 1} but does not surpass 
it. We denote this probability by Q^i(0 | s) and remark 
that, should its dependency with M be known, we would 
immediately be able to discover the desired threshold for 
the number of mutants by adopting s = 1 and specifying 
a lower bound on the probability. In other words, we 
would discover the threshold, call it M*, by specifying 
Q* such that Q^i(0 | 1) > Q* for ah M < M\ This 
follows from the intuitive expectation that Q^i(0 | s) is 
to decrease as M grows for sufficiently large M. 

From its definition as a conditional probability, 
(3^i(0 I s) is given by AB/Qn+i{() \ s), where A is the 
probability that the number of mutants in the network 
eventually increases from s to M and B is the probability 
that, given that it has M mutants, the system eventually 
returns to state without ever increasing its number of 
mutants beyond M. We calculate the values of A and 
B by resorting to discrete-time Markov chains entirely 
analogous to the one we have been using, but now hav- 
ing reduced numbers of states. The first of these chains 
has states 0, 1, ... , M, of which and M are absorb- 
ing while all else remains unchanged. We set A to the 
probability that the system gets absorbed into state M 
having started at state s, that is, A = Pm-\-i{M \ s). 
The second chain has states 0,1,. ..,M + 1, with 
and M -\- 1 the absorbing states and everything else un- 
changed. We set B to the probability that absorption 
occurs at state once the system is started at state M, 
that is, B = Qm+2(0 | M). We then obtain 






i(0|s) = 



Pm+i{M I s)Qm+2(0 I M) 



(13) 



Qn+l(0 I S) 

Closed-form expressions are in general not known for 
Eqs. dU or (fl^ . so we assume that it suffices to consider 
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FIG. 3: (Color online) M as a function of r for Qn+i(0 | 1) 
e. 



the case in which Pv/(lv = ^ for all transient states v. In 
this case, and for r 7^ 1, Eq. (fT3|) yields 



?f+i(0|5) 



(1 - l/r^)(l - 1/r) 1 - 1/r" 

(1 - l/r^)(l - l/r^+iyM \jr - l/r« ' 

(14) 
Readily, Q^i(0 | s) becomes independent of n as n 
grows, regardless of whether r > 1 or r < 1 (but note 
that the two limits differ). For n ^ 00, in Fig. Owe 
show plots of Q^i(0 I 1) = e for different values of e 
and 1 < r < 2. As expected, for fixed e the value of M 
increases with decreasing r in this range, and somewhat 
counterintuitively we see that the rate of increase is ever 
smaller as r approaches 1. However, this is easily con- 
firmed as we realize that, as in the figure, the limit Q of 
Q^i(0 I 1) as r ^ 1 is such that M ^ ^/l/Q. 

Computational results on the distribution of M, the 
maximum number of mutants achieved when the dynam- 
ics ends in extinction, are given in Figs. H] through [6] for 
a variety of topologies. Figure H] refers to the i^-funnel. 
Fig. [5] to random networks with a Poisson-distributed 
number of out-neighbors, and Fig. [6] to a selection of 
networks that includes an instance of each of these two 
types and also the unidirectional ring. The latter fig- 
ure also includes plots of the analytical prediction given 
by Eq. (p!4|) . In this respect, notice that the prediction 
that corresponds to the r value used in the simulations 
matches the data for the unidirectional ring perfectly. 
This, of course, is consistent with the fact that, in this 
case, the requirements of the isothermal theorem are sat- 
isfied and therefore the assumption behind Eq. (p!4|) is ex- 
act. Similarly, in the figure's inset we demonstrate that, 
for the random network in use, the prediction of Eq. (p!4|) 
with a slightly lower value of r is also a perfect match. 

Notice also, in all three figures, that for M = 100 no 
probability is above 10 ~^. This means that, during the 
simulations, of all events that ended in extinction, no 
more than one out of 10^ achieved more than 100 mu- 
tants. As a consequence, should we use M* = 100 (cor- 
responding roughly to Q* = 10~^, given that there were 
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FIG. 4: (Color online) Probability distributions of M for the 
i^- funnel. Data are averages over 10^ events. 




FIG. 5: (Color online) Probability distributions of M for 
random networks with a Poisson-distributed number of out- 
neighbors. Data are averages over 10 events (10 events for 
each of 10^ graphs with at least 0.95n nodes in the GSCC). 
The plot labeled "average" refers to all data for r = 1.1, 
n — 1 100, and z — 4 in all panels. This is compared in 
the same panel to the results corresponding to a single graph 
whose GSCC has 1 045 or 1 070 nodes, for which 10^ events 
have been run. 



at least 10^ events) and count as an event ending in fix- 
ation any event achieving more than M* mutants, we 
would be introducing a deviation of no more than 10 ~^ 
with respect to the actual value of the fixation probabil- 
ity. But the fixation probability obtained by full simula- 
tions of the evolutionary dynamics is itself subject to the 
so-called standard error that is inherent to any Monte 
Carlo simulation. If p denotes the fixation probability 
calculated after N events, then the standard error is the 
standard deviation of the O's (extinctions) and I's (fixa- 
tions) accumulated along the events divided by \/]V, that 
is, y^p{l — p)/N . This function of p is plotted in Fig. [3 
for different values of N ^ along with a flat line for the con- 
stant 10 ~^. Clearly, the additional deviation introduced 



FIG. 6: (Color online) Probability distributions of M for as- 
sorted topologies with r = 1.1. Data are averages over 10^ 
events for the i^-funnel and the unidirectional ring, 10^ events 
for the random network (10^ events for each of 10*^ graphs 
with at least 0.95n nodes in the GSCC). 
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FIG. 7: (Color online) Inherent uncertainty associated with 
the Monte Carlo computations of the fixation probability, as 
given by the standard error of the fixation probability esti- 
mate p. 



by the use of M* = 100 is negligible when compared to 
the standard error. 

Speedup figures resulting from the use of early fixation 
detection for M* = 100 are shown in Figs. [8] and [9l re- 
spectively for the i^-funnel (with several n and K values) 
and the if-layer random networks of [25] (with two values 
of n, K = 5, and several values of the a parameter that 
in [25] is used to control the layer-selection mechanism 
as the network is grown by the addition of new nodes). 
Plots in the latter figure are given against 5'(X, Y"), which 
in [25] is used to indicate, if sufficiently above 0, how 
close each network is to having, like the i^-funnel, an 
exponentially growing number of nodes per layer as one 
moves "upstream" through the layers. Clearly, speedups 
are very significant, particularly for the i^- funnel with 
the largest values of n and K and all random networks of 
10 000 nodes. We also remark that, in Fig. [51 the fact that 
slopes increase with K is closely related to Eq. ^: as the 
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FIG. 8: (Color online) Speedups resulting from the use of the 
M* = 100 threshold on the iC- funnel with r — 1.1. Data are 
averages over 10^ events. 
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FIG. 9: (Color online) Speedups resulting from the use of 
the M* = 100 threshold on the random networks of 25] with 
r = 1.1. Data are averages over 10^ events (10^ events for 
each of 10^ graphs), except for n — 10 000 with a = 4, in 
which case 10 events are used for each of 500 graphs. 



fraction of events that lead to fixation grows with K^ we 
expect from Fig. [T]that the average number of steps be- 
have likewise; qualitatively, this explains why speedups 
increase with K. 



V. FINAL REMARKS 

In Section |T] we mentioned that, in [25|, we were unable 
to extend to larger values of n our conclusions regarding 
the fitness-amplification properties of the i^-layer ran- 
dom networks we used in some of this paper's experi- 
ments. Provisioned with the technique of Section HVl we 
can now bypass the computational difficulties that ham- 
pered our progress in that occasion by employing early 
detection of fixation. Doing this for M* = 100 has re- 
sulted in the data shown in Fig. [TOl from which it is fi- 
nally clear that, also for n = 10 000, it is possible to grow 
layered networks that achieve significant fitness amplifi- 



FIG. 10: (Color online) Results for the random networks of 
^25i] with K — ^ and r = 1.1. Each graph whose layer popula- 
tions correlate with those of the i^-funnel by more than 0.9 in 
the sense of the Pearson correlation coefficient is represented 
by its fixation probability and its S(X^ Y) value. Data are 
averages over 10^ events (10^ events for each of 10^ graphs), 
except for n — 10 000 with a = 4, in which case 10^ events 
are used for each of 500 graphs. Dashed lines indicate p — p\ 
through p3. 



cation. As we see in the figure, for if = 5 many grown 
networks have p values between p2 and pa. We note, 
however, that this is still an easier scenario than that of 
[25|, where if = 10 was used for n = 10 000, since in 
the present case we had to calculate the speedups given 
in Fig. [9] and these required that simulation events be 
carried out to completion. 

Another interesting by-product of our use of the 
threshold number of mutants M* is that, should the dy- 
namics be started with s > M* randomly placed mu- 
tants of equal fitness, then fixation would occur almost 
surely. This is so because the probability of there being 
so many mutants in a dynamics that is bound to extinc- 
tion is as small as allowed by the choice of Q*. However, 
note that Eq. ([H]) is of no immediate help in quantifying 
the "almost surely," since it is conditioned on the dy- 
namics ending in extinction and therefore does not apply 
to those cases in which the number of mutants becomes 
large enough that extinction is unlikely. Nevertheless, 
whenever the isothermal theorem holds, Eq. ([1]) implies 
that fixation from the initial s mutants occurs with prob- 
ability (1 — l/r^)/(l — 1/r"^) (see also [1]), which is asymp- 
totically equal to 1 for r > 1 as both s and n grow. 

It is also worth mentioning that, because this study has 
been targeted at directed graphs, Eqs. ([7|) and (|8|) are also 
applicable to the special case of an undirected graph and 
can lead to useful insight also in this case. Specifically, 
suppose we take any strongly connected D without an- 
tiparallel edges and make it functionally undirected by 
adding to it the antiparallel counterpart of every one of 
its edges. If E and the Pi sets continue to refer to the 
original D^ then the contribution of each mutant / to 
T.(Li)eE ^jit)/\Pi\ in Eq. (J8j) jumps from T.(u)eE Vl^il 



toE 



ii,l)eE 



1/|P,;I 



E 



il,i)eE 



1/|P^|, therefore leading to 



a smaller P^ {t) / P~ (t) ratio. On the other hand, if D 
already has antiparallel edges, then one curious special 
case is that of the 2-superstar of [1], which has a cen- 
tral node with n — 1 peripheral nodes that connect to 
it through antiparallel edge pairs. This graph is already 
functionally undirected and, for n ^ oo, p ^ P2. So 
the 2-superstar is somewhat of an exception with regard 
to the P~^{t)/P~{t) ratio for undirected graphs. In fact, 
thus far we have been unable to find any other undirected 
graphs for which p > p2, but to the best of our knowledge 
the question of whether any exist remains unsettled. 

We note, finally, that fixation for the 2-superstar, when 
it happens, is bound to take a considerable number of 
steps to occur (this can be seen in Fig. (TJ since the 
2-superstar and the 2-funnel with the same number of 
nodes are the same graph). For the 2-superstar, the 
sum from Eq. (|8]) mentioned above is Mp{t)/{n — 1) + 



mc(t)(n — 1), where Mp(t) is the number of peripheral 
mutants at step t and mdt) is either 1 or 0, indicating 
respectively whether or not a central mutant exists at 
that step. Clearly, for r only slightly above 1, obtaining 
P+(t) > P~{t) depends crucially on the existence of the 
central mutant. This, in turn, can easily change from 
step to step until fixation is eventually approached. It 
is precisely in cases such as this that the early estimates 
of the fixation probability introduced in Section |IVl are 
most useful. 
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